Simulating Stochastic Processes
Contents
Simulating Stochastic Processes#
1. Discrete time Markov Chain with finite state space#
Let \(\pi_t\) be the probability distribution at time t where \(t = 0, 1, 2, \dots\). Let \(\mathbb{P}\) be the \(n\times n\) one step transition probability matrix. By the Komolgorov-Chapman theorem, \(\pi_{t+1} = \pi_t\mathbb{P}\) and
\[\pi_t = \pi_0\mathbb{P}^t\]
The limiting probability \(\pi_\infty = \pi_\infty\mathbb{P}\). Thus, \(\pi_\infty\) is the eigen vector of \(\mathbb{P}^{'}\) for eigen value = 1.
Specify the states of the process#
states = c("ATLANTA","CHICAGO","DC","NYC")
nstates = length(states)
Initial Probability at time 0#
pi_0 = c(1,0,0,0)
Transition Probability Matrix#
tranP = matrix(runif(nstates*nstates),nrow=nstates,ncol=nstates)
for(i in 1:nstates) tranP[i,] = tranP[i,]/sum(tranP[i,])
The Probability distribution at time n#
ntime = 20
pi_n = matrix(0,nrow=ntime,ncol=nstates)
tranP_n = diag(nstates)
for(n in 1:ntime){
tranP_n = tranP %*% tranP_n
pi_n[n,] = pi_0 %*% tranP_n
}
print("the probablity distribution at time n")
pi_n
print("the n-step transition probability matrix")
tranP_n
[1] "the probablity distribution at time n"
| 0.4153128 | 0.10983998 | 0.09042471 | 0.3844226 |
| 0.3363779 | 0.06879544 | 0.23445368 | 0.3603730 |
| 0.3286973 | 0.07117348 | 0.24573353 | 0.3543956 |
| 0.3281437 | 0.07119461 | 0.24769770 | 0.3529640 |
| 0.3281463 | 0.07126231 | 0.24781316 | 0.3527783 |
| 0.3281550 | 0.07126842 | 0.24782103 | 0.3527556 |
| 0.3281569 | 0.07126948 | 0.24781987 | 0.3527537 |
| 0.3281572 | 0.07126956 | 0.24781958 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781953 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.24781952 | 0.3527536 |
[1] "the n-step transition probability matrix"
| 0.3281573 | 0.07126956 | 0.2478195 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.2478195 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.2478195 | 0.3527536 |
| 0.3281573 | 0.07126956 | 0.2478195 | 0.3527536 |
Limiting probability distribution#
pi = eigen(t(tranP))$vector[,1]
pi = pi/sum(pi)
pi
- 0.328157274996289+0i
- 0.0712695617263433+0i
- 0.247819519069673+0i
- 0.352753644207694+0i
Simulation#
nsim = 50
x = rep(0,nsim)
x[1] = sample(1:nstates,1,prob = pi_0)
for(i in 2:nsim){
x[i] = sample(1:nstates,1,prob = tranP[x[i-1],])
}
freq = table(states[x])
freq/sum(freq)
plot(1:nsim, x, xlab="time", ylab="state", type="l")
ATLANTA CHICAGO DC NYC
0.32 0.04 0.30 0.34
Animation#
data = rbind(c(1,1),c(2,2),c(2,0),c(3,1))
data = data[x,]
library(plotly)
df <- data.frame(
x = data[,1],
y = data[,2],
f = 1:nsim,
s = states[x]
)
fig <- df %>%
plot_ly(
x = ~x,
y = ~y,
frame = ~f,
type = 'scatter',
mode = 'markers',
color = I("blue"),
size = I(50)
)
t <- list(
family = "sans serif",
size = 18,
color = toRGB("red"))
fig <- fig %>% add_text(x = data[,1],y=data[,2],text = states[x],textfont=t,textposition="top right")
fig <- fig %>% layout(xaxis = list(range = c(0, 4)), yaxis = list(range = c(0,3)), showlegend = FALSE)
fig
Loading required package: ggplot2
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
2. Ransom walks#
A random walk is a Markov chain with infinite state space. Let \(\pi_t\) be the probability distribution at time \(t\). Let be the \(n\times n\) one step transition probability matrix \(\mathbb{P}\). The Kolmogorov-Chapman theorem still holds, i.e.,
\[\pi_{t+1} = \pi_t\mathbb{P}\]
This indicates that \(\pi_{t+1}^i = \pi_t^{i-1}\times q + \pi_t^{i+1}\times p\)
Initial probability distribution#
x[1] = 0
Transition probability#
tranP = function(current,p){
return(ifelse(runif(1)<=p,current-1,current+1))
}
y = 1:10000
for(i in 1:10000) y[i] = tranP(3,0.5)
sum(y==2)
5009
Simulation#
nsim = 1000
t = 100
p = 0.5
x = matrix(0, nsim, t)
for(j in 1:nsim){
for(i in 2:t){
x[j,i] = tranP(x[j,i-1],p)
}
}
hist(x[,100])